Coordinate systems and maps
Goals:
In the first exercise, you use the Grammar of Graphics to fix a map with
tmap()by modifying the provided code. The libraries, data, and initial code is provided below.In the second exercise, you inspect and assess the compatibility of the provided spatial datasets
In the third exercise, you
- extrapolate a missing CRS, and
- define the missing projection in R
- In the fourth (and optional) exercise, practice how to
- reconcile CRS of your spatial data
- reproject raster and vectors to a shared SRS
- In the fifth exercise, improve your a map of Aarhus from Week 2 following the best map design principles.
- plot layers together
- limit your data into an area of interest
- Create a mirror pair of honest and lying maps of Danish 2025 kommunal election results
Exercise 1 - Fix a map of Denmark
In this exercise you will learn to make a map with tmap
library, by adding spatial data layers and modifying the arguments that
specify their rendering
Data sets
We will use two data sets: elevation and
counties. Both originate from the UC Davis global databaset
and you can redownload them with getData(). (see Week02,
Exercise 12) The first one is an elevation raster object for Denmark,
and the second one is a spatial dataframe object with polygons
representing the 99 Danish municipalities.
Existing code
Here is the code to create a new map of Denmark. Your role is to improve this map based on the suggestions below.
# get elevation
elevation <- raster("../data/DNK_msk_alt.grd")
# get counties
counties <- readRDS("../data/gadm36_DNK_2_sp.rds")
# starting a static map in tmap v3
tm_shape(elevation) + tm_raster(title = "elev", style = "cont", palette = "BuGn") +
tm_shape(counties) + tm_borders(col = "red", lwd = 3) + tm_scale_bar(breaks = c(0,
100, 200), text.size = 1) + tm_compass(position = c("LEFT", "center"), type = "rose",
size = 2) + tm_credits(text = "A. Sobotkova, 2024") + tm_layout(main.title = "My terrible map",
bg.color = "orange", inner.margins = c(0, 0, 0, 0))# starting a static map in tmap v4
tm_shape(elevation) + tm_raster(col.legend = tm_legend(title = "elev"), col.scale = tm_scale_continuous(values = "brewer.bu_gn",
midpoint = 0)) + tm_shape(counties) + tm_borders(col = "red", lwd = 3) + tm_scalebar(breaks = c(0,
100, 200), text.size = 1) + tm_compass(position = c("LEFT", "center"), type = "rose",
size = 2) + tm_credits(text = "A. Sobotkova, 2025") + tm_title("My terrible map") +
tm_layout(bg.color = "orange", inner.margins = c(0, 0, 0, 0))Tasks
- Change the map title from “My terrible map” to “Denmark’s municipalities”.
- Update the map credits with your own name and today’s date.
- Change the color palette to yellow to brown sequence (“-RdYlGn” in v3 and equivalent in v4). (You can also try other palettes from http://colorbrewer2.org/ or cols4all::c4a_gui())
- Put the north arrow in the top right corner of the map or some better position.
- Improve the legend title by adding the used units (masl).
- Increase the number of breaks in the scale bar and move it to top right.
- Change the borders’ color of the municipalities to black. Decrease the line width.
- Change the background color to some reasonable color of your choice.
Exercise 2 - Create and inspect spatial data
We will use three data sets: elevation,
places and nitrates, explore them, and make a
map. The first one is an elevation raster object for Denmark you already
used above. The second dataset is a list of places you created in week
1. The third is a geochemical dataset from landbrug.dk.
You have already read in the raster, and so now the focus is on the vector data. You need to create these from googlesheet/csv before you can inspect them. Follow the instructions and answer the questions related.
Preparation: Create vector data from a csv with LatLong and WKT
Load the two datasets:
- read in the places googlesheet using the
read_sheet()function (grab it from W01 exercise) - use
read_csv()to grab nitrates.csv file and useslice()function to grab the first 5000 records upon reading it in to reduce its size.
# Load your DK places googlesheet from W01
gs4_deauth() # if the Gdrive authentication is not working for you
places <- read_sheet("https://docs.google.com/spreadsheets/d/1PlxsPElZML8LZKyXbqdAYeQCDIvDps2McZx1cTVWSzI/edit#gid=1817942479",_______)
# Load the first 5000 rows in .csv-file and save it as "nitrates":
nitrates <- read_csv("../data/nitrate.csv") %>%
slice(1:5000)Convert the tabular data into simple feature using their geometry columns
Tabular data contain spatial data, that you can use to convert these
tibbles into simple features! Note the lat/long columns in
places and wkt column in
nitrates. Inspect these columns and convert the objects to
simple feature using the st_as_sf() function using the
coords and the wkt arguments respectively!
- take
places, filter missing coordinates away and then use thest_as_sf()function with thecoordsargument to cast into a simple featureplaces_sf - for
nitrates, use thest_as_sf()function with thewktargument to convert to simple featurenitrates_sf. - use
plot(places_sf$geometry)andplot(nitrates_sf$WKT)to view each simple feature in space.
Instructions and questions
Type answers to the questions as code comments next to or under the code used
Display the
nitrates_sfandplaces_sfobjects and view their structure. What can you say about the content of each file?- What type of data does it store?
- How many attributes does it contain?
- What is its geometry? (points, lines, polygons, surfaces, network)
- What kind of coordinate system does it use?
Display the
elevationobject and view its structure. What can you say about the content of this file?- What type of data does it store?
- What is the coordinate system used?
- How many attributes does it contain?
- How many dimensions does it have?
- What is the data resolution?
Can you plot the files on top of one another? Why?
Exercise 3 - Define coordinate systems
A machine-readable CRS is prerequisite to spatial visualisation and analysis, but not all spatial datasets that float around have its CRS defined. Here you learn to define/specify the coordinate systems for datasets where the CRS is NA/missing/unrecognized by R. You can fix this situation if YOU know what the CRS should be.
Look at the coordinates in places and
nitratesagain. While both refer to places in Denmark, their
formats are different, signalling different CRS. Can you guess what they
are?
placescontain lat/long columns which you have created on the basis of GoogleMaps which usesWeb Mercator, or WGS84 CRS. This translates into EPSG 4326.nitratescontain WKT column with planar coordinates provided by landbrug.dk server. The source server specifies that all data is ETRS89 UTM32N, which corresponds to EPSG 25832.
Instructions:
- define CRS for each spatial object with
st_set_crs()using the CRS they are compatible with. - write the CRS into the name of these defined object in R and write the CRS into the name!
Your solution and answers
# /Start Code / #
places4326 <- places_sf %>% _______(_________)
nitrates25832 <- nitrates_sf %>% _________(________)
places4326
nitrates25832
# /End Code/ #Wonderful, your vector data now have defined spatial metadata. In the next step, reconcile these two CRSs so you can display them in the same map.
Exercise 4 - Reconce coordinate systems
Now that you know that coordinate systems differ, make them compatible. You have two final CRS choices. Convert your objects to one or another CRS so that you can display the nitrate samples and DK places over the elevation map of Denmark.
Remember that you use st_transform() to reproject vector
data, and projectRaster() on raster data to change their
CRS.
Instructions
Option 1: Transform the
nitratesdataset into the coordinate reference system used in theplacesobject.- Create a new object
nitrateswith theplacescrs. You can label itnitrates_####writing the EPSG out for easy differentiation. - Visualize the results using the
plot()function.
- Create a new object
Option 2: Reproject the
elevationandplacesdata into the coordinate reference system used in thenitratesobject.- Create a new object
elevation_####with thenitratescrs. - Visualize the results
(
elevation_####`` together withnitrates`)
- Create a new object
Exercise 5 - Improve your map of Aarhus from Week 2
Make your Week 2 map publication-ready! Download data about DK
available online (if needed). Make a clean and legible map of Aarhus and
a focus area (see below) with tmap library following the
Exercise no.1 above and the guide here or here
Additional challenge:
- Draw a bounding box (hint
st_make_grid()to generate it) or a buffer around the district you live in or some other area of interest, and generate a more detailed map of this area, while leaving Aarhus + Denmark’s boundary in an inset for orientation. Make sure to add title and an explanation for your choice. Include a legend for any environmental background data. - Choose a reasonable classification if you are visualising environmental or quantitative data of any sort.
- Figure out how you could print this result as .png ( explore
ggsave) so the ratio and position of the inset is preserved
Exercise 6 - Challenge: Make a map that is true and that lies
You are a political consultant. Your client wants a map of Danish election results (2025) that makes their party look strong. Choose your party (e.g. DF) and create:
- An honest map: Follows cartographic best practices
- A misleading map: Technically correct but emphasizes your client’s strength (through color, classification, or geographic focus)
- A 300-word memo: Explains how the misleading map manipulates perception without lying
Get data
# get municipalities
muns <- readRDS("../data/gadm36_DNK_2_sp.rds") %>%
st_as_sf()
# get votes
library(googlesheets4)
gs4_deauth()
votes <- read_sheet("https://docs.google.com/spreadsheets/d/1R77xvSEdEF5UMF6u1POlHVyYHbB1zFkNqEYTkwD7GPw/edit?")
# align municipality names
`%nin%` <- Negate(`%in%`)
votes$Division %in% muns$NAME_2 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE
[25] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[37] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE
[49] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[73] TRUE TRUE TRUE
[ reached getOption("max.print") -- omitted 23 entries ]
[1] "Aarhus" "Nordfyn" "Høje-Taastrup" "Copenhagen"
muns$NAME_2[31] <- "Aarhus"
muns$NAME_2[21] <- "Høje-Taastrup"
muns$NAME_2[25] <- "Copenhagen"
# join up
mun_vote <- muns %>%
select(NAME_2) %>%
left_join(votes, by = c(NAME_2 = "Division"))
mun_voteSimple feature collection with 99 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 8.076389 ymin: 54.55903 xmax: 15.19306 ymax: 57.75153
Geodetic CRS: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
First 10 features:
NAME_2 A B C F I M O V Æ Ø Å Others
1 Albertslund 25.3 11.3 7.4 16.1 3.8 1.2 11.9 3.3 1.5 14.4 2.1 1.8
2 Allerød 14.2 9.5 21.5 12.4 6.1 2.0 8.0 17.7 1.1 5.9 0.8 0.8
3 Ballerup 33.9 7.9 8.9 10.6 6.1 2.4 13.7 5.6 1.6 7.2 1.1 1.1
4 Bornholm 46.0 2.6 5.5 5.6 1.9 0.6 16.2 7.6 2.7 7.1 0.5 3.6
5 Brøndby 28.3 10.2 8.3 10.3 4.2 2.1 17.0 4.2 2.0 10.6 1.6 1.4
geometry
1 MULTIPOLYGON (((12.37129 55...
2 MULTIPOLYGON (((12.20744 55...
3 MULTIPOLYGON (((12.42482 55...
4 MULTIPOLYGON (((14.91719 55...
5 MULTIPOLYGON (((12.44434 55...
[ reached 'max' / getOption("max.print") -- omitted 5 rows ]
Look at the data
Is everything in place, any municipalities missing?
Choose a map that is honest and one that lies
# 'Jenks' style further smooths over the gaps
b <- tm_shape(mun_vote) + tm_polygons(col = "A", style = "jenks", n = 6)
b# quantile style divides into 5 even groups
c <- tm_shape(mun_vote) + tm_polygons(col = "A", style = "quantile", n = 5)
c# Equal interval style divides the distribution into even groups
d <- tm_shape(mun_vote) + tm_polygons(col = "A", style = "equal", n = 5)
d# Write maps above to objects and plot them side by side with tmap_arrange()
# for better comparison
tmap_arrange(a, b, c, d)Refs
Tennekes, Martijn. 2019. Tmap: Thematic Maps. https://CRAN.R-project.org/package=tmap.